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Abstract 

We present a study of the dynamics in — > K~^K~Kg decays with approximately 230 million 
BB events collected by the BABAR detector at SLAC. We find that the Dalitz plot distribution is 
best parameterized with the (pK^ mode, an S-wave K~^K~ resonance near 1500 MeV/c^, and a 
large non-resonant contribution. We set limits on resonances not included in our model, and study 
models for the non-resonant contribution. 
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1 INTRODUCTION 



Charmless three-body decays of the B meson are a rich laboratory for the physics of the Standard 
Model (SM), providing information on both the weak sector and the dynamics of the strong inter- 
action. Decays to the final state K^K^ Kg, which are dominated by 5 ^ sss amplitudes, are of 
particular interest due to their sensitivity to physics beyond the SM. However, little is known about 
the dynamical structure of this decay. Such an understanding is important because although the 
branching fraction is large [1, 2], the final state is a mixture of CP-even and CP-odd states [1, 3], 
hence any measurement of the CP asymmetry is diluted due to the presence of both components. 
A key to more precise measurements of the CP asymmetry is to identify and parameterize various 
contributions to the final state. An amplitude analysis, incorporating interference between decay 
submodes, allows for the development of a model of the decay dynamics and the extraction of the 
partial branching fractions and relative phases for the modeled resonances. 

In this paper we present results from a full amplitude analysis of the decay K'^K''Kg. 
The decay kinematics of a spin-zero particle into three spin-zero daughters are completely de- 
termined by two degrees of freedom. In terms of the invariant masses of daughter pairs with 
four-momenta Pi, mf^ = (pi +Pj)'^, the decay rate is 



1 1 



(27r)3 32M|o 2 



\A\^ + \A\^ 



\A\' 



\A\ 



(Am<ir)2 -I- 1 



(1) 



where M^o, Am^, and r are the mass, mixing frequency, and lifetime of the B^, respectively. A 
(A) is the signal B^ (B^) decay amplitude and gtag is -1 (1) when the other B meson in the event 
IS a Summing over qtag, we obtain 



dr{mK+K-,'mK+Ko) 



(27r)3 32M2o 



1 

2^ 



1^1' + 1^1 



(2) 



Assuming .4. ~ .4, corresponding to no direct CP asymmetry, this expression simplifies to 

C^r(T»A-|A- -'"a-IA-") 1 1 



'A I A»J 



(27r)3 32M|o 



(3) 



This assumption is consistent with the cosine term observed in time-dependent CP asymmetry 
measurements [2, 4], and predicted by preliminary calculations based on the QCD-factorization 
model [5, 6]. 

To describe the decay dynamics, we parameterize the amplitude in the isobar model [7], as a 
sum of contributions A = Xlr Cr- Ar, where is a complex coefficient and the index r runs over the 
resonances in the model plus a non-resonant component. For a resonance r formed in the variable 

2 



m. 



Ar {ml , J = Tr {ml ) x Pb ( ) x P^ ( |pl \) x Zl {m ,Pk), (4) 

where the momenta p arc measured in the resonance rest frame, is the dynamical function of 
the resonance, Pl (Pb) is a Blatt-Weisskopf barrier factor for the resonance {B meson) decay, and 
Zl is a Zemach tensor [8, 9]. 
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The Blatt-Weisskopf factors Fl{z), with z = \pi\R and the range R = 1.5 GeV ^ [10], depend 
on the resonance angular momentum L [8] : 




(5) 

_ I g + Szj + zl 

^^=^^^) - V9 + 3Z2 + .4- 

Here zq = z{\pi\Q), where {pHq = \pi\ when rriij = m^, the resonance mass. For the factor Fb for 
the parent B meson, we use a range R of zero. The angular distribution of the decay products for 
different L is given in the Zemach tensor formalism by [9] : 

ZL=o{Pi,Pk) = 1, 

ZL=i{Pi,Pk) = -^Pi-Pk, (6) 
16 r 

ZL=2{Pi,Pk) = y [HPi-Pkf - {\Pi\\Pk\] 
Most resonances are parameterized with the relativistic Breit-Wigner form 

Triml) = — 2 ^ w 

— nij^j — irriri. {mij) 

r(mij) is the mass-dependent width, given by 

where is the nominal resonance width. For the /o(980), we use the coupled-channel parameter- 
ization of Flatte [11], with couplings taken from recent measurements [12, 13]. 

The inclusive charmless branching fraction B{B^ K'^K~K^) is measured to be (24.7ib2.3) x 
10~^ [1, 2], where (4.2^° 5) x lO^*' is from the cpK'^ channel [14, 15, 16]. We have previously reported 
on an angular moment analysis for this decay [3]. As a function of rrij^+x-, we found the (p to be 
the only statistically significant P-wave component. The S-wave contribution is spread across the 
kinematic range, with a peaking structure at around 1500 MeV/c^. Contributions of higher waves 
were consistent with zero. Our present analysis model is motivated by these results. In addition 
to the <^(1020)i^° channel we include an S-wave resonance with a mass of about 1500 McV/c , 
the X(1500). An additional S-wave component, the "non-resonant" (NR) contribution, is added 
to account for the S-wave events spread across the phase space. We also try adding the /o(980) to 
the fit model. 

The absence of higher waves in the mx+x- angular analysis suggests the lack of resonant 
channels decaying to K~^Kg. The Oq (980)i^''" branching fraction has been measured as 

small [17]. Additionally, there are theoretical reasons to doubt the possibility of significant isovector 
contributions [18, 19]. 

Several decays involving b ^ c transitions can contribute to the K^K^K^ final state. We 
coherently add an amplitude for the XcoKg mode to our isobar model, while the and Df 

K" channels are added incoherently and are described with Gaussian shapes to account for the 
finite resolution in m{K'^K~) mass. 
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Unfortunately, the present statistics do not allow for the parameterization of the P-wave outside 
of the (j)Kg region and the subsequent extraction of the relative S — P phase. Hence, we focus on 
the parameterization of the S-wave contribution which dominates the decay [3] . Theoretical models 

of the NR amplitude [5, 20] do not reproduce the distribution observed in data. Without reliable 
theoretical guidance, we test several ad-hoc parameterizations. The three most successful models 
are: 



/ 2 

rriK+K- log ' 



(9) 



1 + /3„,e">«m^+;^,, and (10) 
exp [pnrm^K^K^ , (H) 

multiplied with isobar amplitudes and constant phases, and the parameters f3nr are free in the 
fits. The first parameterization is inspired by a model from Ref. [5], where to improve agreement 
with data we have dropped terms involving b u transitions and those that give a dependence on 
m^+^o. The second parameterization is a linearization of theoretical models [20], and the third 
model was successfully used in an analysis of the decay S+ — K''^K'^K~ [21]. We choose the 
model of Eq. (11) to parameterize the NR amplitude. This model gives good agreement with data 
and allows for direct comparison of the Pnr parameter with that determined in — > K~^K~^K~ 
decays. 

We also test for a possible NR amplitude dependence on the K^Kg mass by adding terms with 
a linear dependence on m\+K° ■ 

Therefore, our main model includes four coherent contributions: the (/)(1020) {m = 1019.456 MeV/ c^, 
r = 4.26MeV/c2), X(1500) (m = ISOTMeV/c^, T = lOOMeV/c^), and XcO {m = 3415.19 MeV/ c^, 
r = 10.1 MeV/ c^) decaying to K+K'- and the NR S-wave [22]. 



2 ANALYSIS METHOD 

This analysis uses a data sample of approximately 230 million BB events collected with the BABAR 
detector [23] at the SLAC PEP-II e+e^ storage rings operating at the T{AS) resonance. The basic 
method used for reconstruction and selection of S°-candidates is described in Ref. [2]. We charac- 
terize events with standard topological variables which distinguish between the jet-like structure of 
the dominant continuum e'^e~ qq {q = u,d. s,c) background and relatively isotropic B decays. 
These variables are combined into a Fisher discriminant J- [2] . The shape of the J- distribution for 
continuum events is seen to vary as a function of the Dalitz plot location, becoming more "signal- 
like" away from the Dalitz plot edges. As a result, we do not include it in the fit, and instead 
impose a requirement on the value of to select events for the fit. 

Two kinematic variables are used to characterize B^ candidates, the energy difference AE = 
Eb — \/s/2 and the beam-energy-substituted mass ttt-es = y^(s/2 + pi -pbY/E^ — p^. Eb is the 
reconstructed energy of the B^ candidate in the CM frame and i/i is the total CM energy, (p,, Ei) 
is the initial e^e~ four-momentum and pB is the reconstructed S'^-candidate momentum, both 
measured in the laboratory frame. For signal events, /S.E peaks at zero with a resolution of 18 MeV, 
while rriES peaks at the B^ mass with 2.6 MeV/ resolution. We initially retain candidates with 
lAE'l < 200 MeV and mss > 5.2GeV/c^. For the Dalitz plot fit, we keep candidates in a signal 
region (SR) of lAE'l < 60 MeV and mss > 5.26GeV/c^. We also define a sideband (SB) region, with 
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|A£^| < 200 MeV and 5.2 < uies < 5.26GeV/c^, to study continuum backgrounds. For calculation 
of the Dalitz plot variables, candidates are refit with their mass fixed to the world average value [22], 
constraining the B° candidate to the kinematically-allowed Dalitz plot region. 

Backgrounds due to B decays are studied with samples of simulated events. The largest con- 
tribution arises from the random combination of tracks from both B mesons. This combinatorial 
background is difficult to differentiate from continuum background, and is accounted for in the fits 
with the continuum background. B decays with a missing pion in the final state combination fall 
outside both the SR and SB. 

We use a high-statistics sample of simulated signal events to evaluate the signal efficiency in bins 
on the Dalitz plot. The efficiency varies primarily as a function of mx+K-i falling from roughly 
35% for low rux+K- to about 10% at high mx+K-- The average efficiency is 26%. Signal events 
which arc misrcconstructcd with a daughter from the other B meson give a negligible contribution, 
making up less than 0.5% of the total events. 

An unbinned extended maximum likelihood (ML) fit to the events in the SR is used to extract 
event yields and amplitude coefficients. Parameters of the probability density functions (PDF) for 
signal AE and mss are determined in simulated events, and those for the continuum background 
are fit to data. For the amplitude fit, we parameterize the event kinematics with the variables 
rriK+K- and cos 6h, where Oh is the angle between the fiight direction of the and Kg in the 
K'^K~ center-of-mass (CM) frame. This change of variables introduces a Jacobian term \J\ = 
{2mp^+ x-){'^\'Pi\\Pk\) in the signal PDF, where {\pk\) is the (Kg) momentum in the K'^ K~ 
CM frame. The total PDF is formed as V{mEs) • V{AE) ■ ^'(Dalitz), where for signal events 

P(Dalitz) = c?r(mii-+x-,cos^if) |J| £{mK+K- ^^osOh) (12) 

before normalization, and e is the efficiency. The continuum ^-'(Dalitz) PDF is modeled with a 
histogram of events from the SB region. Bin sizes are variable across the Dalitz plot to account for 
both narrow features and sparse statistics in different regions. 

To explore the multiple solutions possible in the Dalitz fit, we performed several hundred fits 
with the initial amplitudes and phases randomized within reasonable values. In the fits, the am- 
plitude and phase of the NR contribution are fixed. The relative contributions of each component 
are evaluated with the fit fraction FFj., defined ignoring the interference between amplitudes, as 

pp ^ / \CrAr?\J\dmK+K-dcOSeH 

J\Y,s(^s-As\'^\J\dmK+K-dcoseH' 

3 SYSTEMATIC STUDIES 

We study several potential sources of systematic uncertainty. Fitting to events generated with the 
PDF and with a full MC sample, we test for any fit bias. No significant biases arc found, and we 
conservatively include the largest deviation from the generated parameters as systematic errors. 
The parameters of the PDFs for the kinematic variables are fixed in the fit. We vary them by 
one standard deviation and take the change from the nominal fit as the systematic error. The 
masses and widths of all resonances except the /o(980) arc taken to be the world averages [22], and 
we vary the values by their listed error and assign systematic errors based on the changes in the 
fit results. For the /o(980), the systematic error due to the parameterization is derived from the 
spread between the solutions given by the coupling values measured by BES and E791 [12, 13]. We 
derive an analogous uncertainty for the X(1500) by fitting with the parameters given by Belle for 
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the /x(1500) — > K~^K~ observed in — > K~ decays [21] and comparing the results with 

the nominal fit that uses the world average values [22]. We estimate the systematic uncertainty due 
to the parameterization of the NR component from the difference in results for the three models. 

An additional contribution to the systematic error arises from the spread of values at each of the 
quoted solutions. We assign the RMS of this spread as the uncertainty. We observe that several of 
the fit parameters are the same for multiple clusters of solutions, with only one or two parameters 
differing. In these cases, we quote the average as our result for the parameters that are nearly the 
same and take the difference from the average as a systematic error on the result. 



4 RESULTS 

There are eleven free parameters in the fit: signal and background yields, the NR parameter jdnr, 
and five strengths and three relative phases. The phases for the three interfering components are 
defined with respect to the NR component for which the phase of the coupling to the B meson is set 
to zero. In the extended ML fit with a total of 1842 events, we find 530 it 28 signal events. Figure 1 
shows the data overlayed on mss ^-nd IS.E projections of the fit function. Signal and background 
yields are correlated with the Dalitz plot parameters at the 1% level. 




Figure 1: Projections of the fit function for mss (left) and IS.E (right) shown with data (points). 
The solid curve is the total PDF, and the continuum background PDF is shown as a dashed curve. 

In several hundred fits with random initial values for parameters varied in the fit, we find two 
two-fold ambiguities resulting in four solutions with similar likelihood values. We list these solutions 
in Table 1, where the fit fraction FFr for each component r was defined in Eq. (13). Note that 
interference terms are neglected, so the sum of fit fractions does not necessarily add to 100%. Plots 
showing projections of the fit function with data are shown in Figure 2. 

We observe ambiguity in the fraction and the phase of X(1500)Xg decays. The solution with 
the relative phase close to zero has a small value for the X(1500)i^g fraction, but the fraction is 7 
times larger when the phase is close to 7r/2. In Tab. 1 and Fig. 2, we denote the former solution 
"A" and the latter "B". Another ambiguity occurs in the 4>K^ component where there are two 
solutions approximately 3.5 radians apart in phase, but consistent in fraction. Here we use "1" and 
"2" to label the solutions. Both ambiguities are reproduced in Monte Carlo studies where a fit to 
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a sample of events generated according to one of the solutions leads to the same ambiguities in the 
fitted parameters as observed in the fit to data. 

All NR models given by Eqs. (9)-(ll) result in a good fit, with consistent fractions and phases. 
Using Eq. (11), we get Pnr = — 0.14ib0.02, which is consistent with the parameter from the linearized 
model of Eq. (10). In the alternative model described with Eq. (9) we fix Pnr = 0.3 GeV/c^ as 
suggested by [5], and we get = 0.5 it 0.1 GeV/c^ if varied in the fit. All fit solutions show the NR 
parameter to be inconsistent with a "flat" (phase-space) model. When we attempt to fit with the 
flat model, the fit does not converge in most cases (99%), and we do not observe distinct clusters 
of solutions. We also probe for a NR dependence on the K^Kg mass by introducing an explicit 
dependence on K^Kg mass into our linear model, but the corresponding coefficients are consistent 
with zero. 

As a goodncss-of-fit measure, we divide the phase space into bins and compute the difference 
between signal-weighted events [24] and the expectation based on the fitted signal model. All bins 
have at least 10 events. We find for a two-dimensional binning the value of xVdof = 77.1/70, and 
in projections in rnx+K- ^"^^ cosOh, respectively, we get x^/dof = 48.9/40 and x^/dof = 11.7/16. 
We find negligible differences in x^/dof among solutions. 

Table 1: Dalitz plot fit results. Labels used in the final column are described in the text. The 
upper limit for the /o(980)Kg mode is at 90% confidence. 



Mode 




Fit Fraction FFr (%) 


Phase 


Solutions 


(^(1020) 


0.085 ± 0.009 


15.4 ±3.4 ±0.6 


-1.47 ±0.16 ±0.2 
2.06 ±0.18 ±0.2 


lA, IB 
2A, 2B 


X(1500) K'^ 


0.71 ± 0.15 


5.2 ±2.2 ±0.9 


-0.17 ±0.21 ±0.1 


lA, 2A 


1.93 ±0.18 


38.9 ±7.3 ±0.9 


1.19 ±0.08 ±0.1 


IB, 2B 


Non-resonant 


7.07 (fixed) 


70.7 ± 3.8 ± 1.7 


(fixed) 


All 


XcoKl 


0.32 ±0.07 


3.1 ± 1.6 ±0.8 


0.2 ±0.4 ±0.6 


All 


/o(980) 


1.57 ±0.41 


5.7 ± 3.2 ± 1.0 (< 9.7) 







5 DISCUSSION 

We find that the majority of the decays to K~^K~Kg belong to the K~^K~ S-wscve that we 
identify as the non-resonant contribution. We find three equally good paramcterizations that 
describe this contribution as an S'-wave with a dependence on K^K^ mass. When we apply the 
same model as in K^K^K~ decays [21], we find consistent values for the shape parameter. 

We do not observe any variation of this component with K^K^ mass. This is consistent with our 
previous angular moment analysis [3] that did not find higher angular moments that would arise 
from such a variation. We find good agreement with a theoretical model [5] only when ignoring 
amplitude terms with features that are not found in the data. 

In addition to the non-resonant S'-wave, we identify a contribution with a K^K~ mass around 
1500 MeV/c^. A two-fold ambiguity is observed in the fraction, depending on the relative phase 
with respect to the large non-resonant background. The nature of this contribution is unclear. 
Identification of this state as the /o(1500) is inconsistent with the measurement of the B^ 
/o(1500)ii'°,/o ^ TT+TT- decay [25]. Since the ratio r(/o(1500) ^ K+K-)/r{fo{1500) tt+tt") 
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Figure 2: Top row: rrix+K- (left) and rrij^+j^o (right) projections of the nominal solutions (lines) 
with data (points). For the fn^+j^o plot, there are two entries per event. The peak at m^+j^o ~ 

1.8GeV/c^ is due to D^K~ and DfK~ decays, while the peaks at higher values of ruj^+j^o are 
reflections of forward and backward 4> decays. Bottom row: Left, cos Oh projection of the nominal 
solutions (lines) with data (points); right, Dalitz plot showing distribution of signal events. Lines 
show solutions lA (black solid), IB (red dashed), 2 A (green dotted), 2B (blue dash-dot) (See 
Tab.l). All plots show background-subtracted signal events [24]. 



is 0.1845 ± 0.0195 [22], we would expect only about two events in B° X(1500)K°. Our dataset 
has insufficient statistics for an independent determination of the mass and width of the state, so 
we use the world average values for the /o(1500) [22]. 

In the low K+K' mass region, we find the P-wave contribution from (pK^ decays with a fraction 
consistent with previous measurements [22]. With the available statistics, we do not observe a 
significant 5- wave under the (/> resonance. We perform a fit with the /o(980) resonance included in 
the model and parameterized with the Flatte lineshape, and set an upper limit on the fraction for 
this decay of 9.7% at the 90% confidence level. 
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6 SUMMARY 



We present the first study of the decay dynamics in K^K~Kg decays and report prehminary 

results on the fractions and relative phases of intermediate states that contribute to the decay. 
This is of particular importance for the interpretation of measurements of CP asymmetries in these 
decays, and searches for physics beyond the Standard Model. We assume no direct CP asymmetry, 
allowing us to use all events in this statistics-challenged analysis. This assumption is supported by a 
cosine term in time-dependent CP asymmetry measurements consistent with zero, and preliminary 
calculations based on the QCD-factorization model [5, 6]. 

One of the main goals of this paper has been to provide input to theoretical studies that can 
result in better models. We extended our previous analysis [3] by providing a parameterization for 
the dominant S-wave contribution, and reported on features that are not seen in our dataset. 
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